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Chapter  1:  Introduction 


The  reduction  of  size  and  cost  for  computing  and  electromechanical  components  is  driving 
the  deployment  of  networked  cyber-physical  systems  (CPS).  Defense,  healthcare,  infrastructure, 
and  transportation  industries  are  building  CPS  to  create  new  services  for  their  consumers.  As  these 
technologies  proliferate  and  are  composed  to  create  vast  systems  of  systems,  there  will  be  an 
increased  risk  of  instabilities  and  cascading  failures  due  to  improperly  designed  control  modules. 

Some  recent  examples  of  control  induced  failures  include  power  grid  blackouts  and 
shutdowns  [7,15]  due  to  faulty  actuators  and  human-interaction.  Large-scale  cloud  applications 
also  suffer  outages  because  of  improperly  designed  software  control  logic  [2],  Even  simple,  stable 
automatic  controllers  may  cause  issues  when  chained  together  with  other  systems,  e.g.  [4,  18,  21, 
23].  These  examples  indicate  that  there  is  a  fundamental  issue  with  composition  at  all  levels  of 
system  abstraction. 

For  the  past  thirty  years,  there  has  been  considerable  effort  to  model,  compose,  and  verify 
pure  software  systems  such  that  they  behave  as  specified  by  a  customer  or  designer.  One  way  to 
deal  with  software  system  composition  is  to  apply  software  modeling  and  analysis  techniques  that 
ensure  the  software  operates  according  to  a  specification  [1],  This  model  checking  process  requires 
a  formal  model  for  each  software  component,  typically  some  finite  state  machine,  and  a 
specification  language  (e.g.  LTL  [19]  or  CTL  [5])  that  captures  system  requirements.  An  analysis 
tool  consumes  these  two  components  and  applies  an  algorithm  to  determine  if  the  desired 
specification  was  met  by  the  software  component.  More  sophisticated  models  of  software 
components  were  developed  to  handle  concurrent,  communicating  systems,  such  as  the 
communicating  sequential  process  language  [10]  and  calculus  of  communicating  systems  [14], 
These  new  modeling  paradigms  introduced  the  notion  of  bisimulation  to  determine  the  semantic 
equivalence  of  one  process  to  another  process. 

Since  CPS  have  major  software  components,  they  inherit  the  same  model  checking 
challenges  as  pure  software,  but  make  the  problem  more  complex  through  interactions  with  sensor 
and  actuator  hardware  as  well  as  heterogeneous  communication  systems.  More  recent  research  has 
built  new  CPS  software  analysis  and  synthesis  tools  that  target  individual  components  of  the 
system.  One  important  extension  from  the  model  checking  community  is  the  development  of 
approximate  bisimulation,  e.g.  [8,  9,  24],  This  work  allows  the  comparison  between  dynamical 
systems  to  determine  if  they  have  similar  behavior  within  some  error  bound.  Approximate 
bisimulation  has  been  applied  to  software  controller  analysis  [11]  and  synthesis  [13]. 

The  contribution  of  this  paper  is  the  TeamBlocks  framework  that  facilitates  the  correct 
construction  of  cyber-physical  systems  from  source  code  to  collections  of  systems.  We  model 
software  and  systems  using  a  modified  hybrid  automata  and  leverage  recent  advances  in  interface 
controllers,  e.g.  [6],  to  link  together  modules  at  each  step  in  the  abstraction  hierarchy.  Furthermore, 
approximate  bisimulation  is  used  to  generate  behavioral  certificates  for  implemented  controllers 
and  compositions  across  abstraction  layers  such  that  they  approximately  match  their  theoretical 
design.  By  integrating  these  theoretical  tools  across  each  abstraction  boundary  in  the  CPS 
hierarchy,  TeamBlocks  guarantees  the  correct  operation  of  the  composed  system  within  user 
specified  bounds,  or  it  detects  that  the  composition  fails. 
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Section  3  introduces  TeamBlocks’s  core  concepts  and  models.  Section  4  discusses  our 
modeling  assumptions  and  applies  the  work  in  Section  3  to  model  software  and  systems.  The 
construction  of  behavioral  certificates  is  presented  in  Section  5  and  the  simulation  results  for  a 
candidate  pair  of  systems  are  demonstrated  in  Section  6.3.  We  conclude  with  a  discussion  of  our 
simulation  results  and  potential  research  extensions  in  Section  9. 
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Chapter  2:  Use  Cases 


Before  proceeding  with  formal  definitions,  it  will  be  helpful  to  point  out  a  few  of  the  use  cases  that 
we  have  in  mind  for  TeamBlocks ’  software  tools,  and  for  the  development  of  the  underlying  theory, 
so  that  we  can  see  where  we  are  headed. 

Example:  Generate  model  from  C++  source  code.  Suppose  we  have  developed  a  controller  to 
be  implemented  on  an  embedded  system.  The  implementation  has  been  written  in  C++.  Now  we 
would  like  to  perform  further  analysis  of  the  system,  based  on  this  C++  controller  implementation, 
including  whatever  quirks  it  may  have.  To  achieve  this,  our  first  step  is  to  generate  a  model  of  the 
C++  implementation  that  can  be  used  with  the  rest  of  the  analysis  tools.  For  a  supported  subset  of 
C++  code,  TeamBlocks  provides  exactly  this  automatic  model  extraction  capability.  Give  it  code; 
it  returns  a  hybrid  automaton  model. 

Example:  Compose  and  execute  models.  Suppose  next  that  we  are  interested  in  assembling  larger 
systems  out  of  smaller  ones,  so  that  we  can  understand  how  the  larger  system  (or  system  of 
systems)  behaves.  TeamBlocks  provides  tools  to  do  this  as  well,  and  to  prove  properties  of  that 
larger  system.  That  is,  TeamBlocks  provides  tools  for  composition  and  execution  of  hybrid 
automaton  models  formed  from  smaller  ones. 

Example:  Check  execution  against  specification.  If  we  have  a  specification  for  how  a  system 
ought  to  behave,  we  would  like  to  test  whether  a  real  system  actually  satisfies  that  specification. 
TeamBlocks  provides  a  Runtime  Validator  that  does  this  for  a  class  of  specification  described  using 
Linear  Temporal  Logic  (LTL). 

Example:  Use  provided  models.  If  we  simply  wish  to  model  a  system-of-systems  without  starting 
from  scratch,  TeamBlocks  provides  a  small  library  of  useful  models  for  example,  of  vehicles  that 
can  be  used  within  the  framework.  These  are  described  using  the  same,  unified  representation  as 
software  and  other  systems. 

Example:  Verify  abstraction.  As  systems  are  composed,  the  dimensions  of  their  state  spaces  are 
multiplied,  and,  as  a  result,  one  often  runs  quickly  into  problems  of  extremely  large  state  spaces. 
We  use  abstraction  to  combat  this  complexity,  by  verifying  that  a  complicated  system  actually 
behaves  like  another,  simpler  system,  to  within  an  approximation  bound.  Given  two  systems,  we 
would  like  an  automatic  proof  that  the  more  complicated  system  actually  does  behave 
approximately  like  the  simpler  one  from  an  input-output  perspective. 

As  part  of  TeamBlocks,  we  have  investigated  problem  formulations  and  developed  prototypes  that 
automatically  search  for  proofs  that  the  one  system  does  in  fact  approximate  the  other. 

With  these  motivating  examples  in  mind,  we  in  the  next  chapter  jump  into  a  description  of 
TeamBlocks ’  underlying  modeling  framework. 


Approved  for  Public  Release;  Distribution  Unlimited. 
3 


Chapter  3:  Modeling  Framework 
3.1  Generalized  Hybrid  Automata 

TeamBlocks  requires  a  modeling  paradigm  that  can  describe  both  control  software,  which  typically 
has  discrete  time  and  state,  and  physical  systems  (e.g.,  the  aerodynamics  of  a  UAV),  which  are 
typically  modeled  with  continuous  time  and  state.  Our  approach  is  to  use  a  single  hybrid- 
automaton  representation  for  both.  In  particular,  we  use  a  kind  of  controlled  hybrid  automaton  with 
output  map,  which  we  refer  to  throughout  as  the  generalized  hybrid  automaton  (GHA).  Its 
definition  follows: 

Definition  1.  Let  a  generalized  hybrid  automaton  (GHA)  be  defined  as  the  tuple 

A  =  (Q,X,U,YXV,H,E,®,R) 


consisting  of: 

•  a  finite  set  Q  =  {ql,...,q\Q\}  of  modes, 

•  a  set  X  =  {Xq}q<=Q,  where  Xq  C  RJ</,  dq  e  Z+ 

•  a  set  U  C  Rm  of  admissible  control  inputs, 

•  a  set  Y  C  Rp  of  possible  outputs, 

•  a  finite  set  E  of  synchronization  labels, 

•  a  set  F  =  {fq}qsQ of  vector  fields  fq :  Xq  x  U  — »  Xq, 

•  a  set  H  =  {hq}qsQ  of  output  maps  hq :  Xq  — >  Y , 

•  a  set  E  C  Qi  x  Q  x  2  of  labeled  edges  where  each  edge  e  =  (s,d,a)  e  E  has  the  following 
components  , 

•  a  set  O  =  {(pe}esEoi guard  functions,  (pe :  Xs x  U  — >  Rdim<I>‘/Xs' U) 

•  a  set  R  =  {Re}  e<BE  of  reset  functions,  Re :  Xs  x  U  — *  Xd 

•  an  initial  condition  Wo  C  { (q,x)\q  e  Qtx 

The  GHA  state  space  is  the  set  W  =  {{q,x)\q  e  Q,x  ^Xq).  We  say  that  an  edge  e  ^ £  is  active 
at  time  t  whenever  all  elements  of  <pe(x(t))  >  0  during  the  execution;  the  system  may 
nondeterministically  transition,  via  its  reset  map,  along  any  edge  whenever  that  edge  is  active. 


1  We  use  the  subscripts  s  and  d  as  shorthand  for  source  and  destination  modes. 
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An  execution  of  a  GHA  is  a  function  w  =  (q,  x):  M+  W,  such  that  w(0)  G  W0,  and  there 
exists  a  sequence  of  times  T  =  {t1(  t2, ... }  <=  K+  such  that  for  each  interval  [tk,  tk+1)  G 
{[0,  ti),  [t±,  t2 ), ... }  the  following  equations  hold  for  all  t  G  [tk,  tk+1): 

ftk+i 

x(t)  =x(tk)+\  /(q(T))(x(T),u(r))dT 

Jtk 

q  (0  =  q(tk) 

0(s  d)(x(t),  u(t))  <  0,  V(s,  d)  G  E  s.  t.  s  =  q(t) 

For  each  tj  G  T,  there  exists  (s,  d)  G  E  such  that  q(tk )  =  s,  q(tfc)  =  d,  x(tk)  = 

R(s,d)(x(.tk),u(tk)),  and  4>(s,d)(.x(tk),u(tk))  =  0.2 

We  define  polynomial  GHA,  denoted  Ap,  as  the  subclass  of  GHAs  for  which  the  mappings 
within  F, H,  O,  and  R  are  finite  polynomials,  and  Wo  is  a  sub-level  set  of  a  finite  polynomial. 

3.2  Composition  Operators 

Asynchronous  Composition  of  GHA 

The  TeamBlocks  analysis  operations  require  the  definition  of  an  asynchronous  product  for  two 
given  GHA.  We  adapt  the  standard  asynchronous  product  between  two  discrete  automata  (e.g.  [3, 
17]): 

Definition  2.  Given  two  GHAs, 

Ai  =  (&,  Xlt  Ult  Ylt  T1;  F1(  H1(  Elt  <J>lf  Rlf  Wli0) 

^2  =  (Q.2>  ^2'  X2,  F2,  H2,  £"2,  @2’  H^2,o) 

their  asynchronous  product,  At\  \A2,  is  composed  of 

•  mode  set  Q  =  x  Q2 

•  continuous  state  space  X  =  {Xlqi  ©  X2,q2)(q1,q2)eQ  where  is  ©  the  direct  sum. 

•  labels  X  =  X1  U  X2 

•  initial  condition  W0  =  {(q,x)\q  G  Q10  x  Q2  0,x  G  X10  ©  X2  0 } 

•  edges  from  the  union  of  three  sets: 

E12  =  {((Si,  s2),  (d  1,  s2),  o )  I  (slt  d±>  a)  EE1A<j  g  2^  YT2  a  s2  e  Q2}  (Advance  Ax 
but  not  A2) 


2  For  notational  convenience,  we  let  g(tk  )  denote  limT^0,g(t  —  r)  for  any  given  function  g. 
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%2  =  {((S1<  ^2)'  (Si,  d2),  o)  I  (s2,  d2,  a)  E  E2  Acr  E  TAA  A  sx  G  Q x}  (Advance  i42 
but  not  A  A 

£”12  =  {((%»  s2),  (di,  d2),  ff)  I  (■%,  dlf  o)EE1  A  (s2,  d2,  a)  6  F2)  (Advance  ^  and  42 
together) 

guard  functions  for  the  three  cases, 

-  A1  guard  condition: 

0e((*i.*2))  =  0i,(s1,d1,<x)O 1)  Ve  £  Fj2 

-  j42  guard  condition: 

0e(C'T>'*-2))  =  02,(s2,d2,<7)  (A2)  ^  ^12 

-  Conjunction  of  v42  guard  conditions: 

~  0e(C*-l<  ^2))  =  (01,(s1,di,<7)  (Al)<  02,(s2,d2,cr)  (^2))^®  ^  ^12 

and  reset  maps  for  the  three  cases, 


-  f?e  ((xlf  x2))  =  (i?lj<7l  (Xi),  x2)  Veefj 2 


-  i?e  ((xlf  x2))  =  (xlf  R2M2  (x2))  Ve  G  %2 

-  fleCCXi.Xz))  =  (Ril<7l(xi)’  R2,q2  (*2))  Ve  G  £12. 

When  a  pair  of  GHA,  v4i,/12,  have  inputs  and  outputs,  say  U1,  U2  and  Ylt  Y2  respectively,  we 
need  to  take  care  that  the  asynchronous  composition  properly  handles  the  dimensionality 
requirements  connecting  these  systems.  To  handle  this  operation,  we  develop  a  feedback 
operator  from  one  GHA  to  another. 

Feedback  Operator 

Given  a  product  of  automata,  we  will  often  want  to  define  the  connection  of  its  outputs  to  its  inputs. 
For  this  purpose,  it  will  be  useful  to  define  a  class  of  simple  selection  operators  H{a, :  R;  — ►  R/  , 
k  <  j,  that  return  a  subset  of  a  vector’s  elements.  With  this  operator,  we  define  a  (partial)  feedback 
map  to  be  a  function  Ki,j  :  YhY  —*■  UjU,  I/I  =  I/I,  that  maps  components  (specified  by  index  set  I)  of 
the  output  in  Y  ,  to  corresponding  components  (specified  by  J)  of  the  input  in  U.  Given  an 
automaton  A  and  a  feedback  map  Kij,  we  get  new  vector  fields  at  all  modes  by  composing  fq  and 
K  as  well  as  new  passports  and  transforms  by  composing  <pe  and  Re  with  K. 


3.3  Timing 

The  hybrid  automaton  model  defined  in  the  previous  section  can  describe  a  variety  of  event- 
triggered  systems.  As  a  special  case,  this  includes  discrete-time  digital  control  systems  with  regular 
sample  intervals.  In  this  section,  we  describe  how  this  is  modeled  by  means  of  a  clock  automaton 
(Figure  1).  Clocks  and  time  will  appear  in  more  detail  in  the  next  chapter. 
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0-1 

4>l (t,  x)  =  x  —  1 
R\(t,  x)  =  ( t ,  x) 


fl(t,  x)  =  (1, 2p)  f2(t,  x)  =  (1,  -2 p) 


a2 

<fo{t,x)  =  -X 
R2(t,x)  =  (t,x) 


Figure  1  The  clock  automaton  is  a  polynomial  hybrid  automaton  that,  when  composed  with  a 
controller  model  via  asynchronous  product,  defines  the  passage  of  time  and  the  sample  rate  p.  Its 
modes  {1,2}  =  Q  have  continuous  state  space  X1=X2=X  =  txt9(t,r),  where  t  is  the 
physical  time,  and  x  is  the  state  of  a  triangle  oscillator. 
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Chapter  4:  Modeling  Software  Dynamics 

TeamBlocks  uses  the  polynomial  variant  of  Definition  1  to  model  control  software.  This  approach 
is  similar  in  essence  to  the  work  of  [22];  however,  by  using  the  full  generality  of  hybrid  automata, 
it  is  possible  to  leverage  approximate  bisimulation  for  validation.  This  section  describes  and 
illustrates  how  polynomial  GHA  model  software  components  as  well  as  how  the  models  are 
automatically  generated  from  C/C++  source  code. 

4.1  Modeling  Assumptions 

Figure  2  presents  the  assumed  control  architecture  that  runs  software  controllers.  It  is  assumed  that 
software  controllers  are  non-blocking  functions  deployed  within  a  timed  process  that  executes 
every  T  seconds.  However,  operations  that  access  hardware,  such  as  the  A/D  or  D/A,  may  be  a 
blocking  operation  through  their  software  APIs.  The  sample  period  sets  the  rate  at  which  the 
system  A/Ds  digitize  the  output  signal  for  control  code  consumption.  Once  the  newly  sampled  data 
is  provided  to  the  software  controller,  it  executes  as  quickly  as  possible  according  to  the  system’s 
processor  clock.  The  control  algorithm  computes  a  new  value  that  is  written  to  the  D/A  for 
conversion  into  an  actuator  signal. 


Figure  2  A  diagram  of  the  assumed  control  software  execution  when  deployed  on  the 

cyberphysical  system. 

Modern  software  controllers  are  typically  implemented  on  embedded  systems  with  short 
instruction  periods  (jjs  to  ns)  using  a  zero-order-hold  strategy.  For  example,  an  unmanned  aerial 
vehicle  autopilot  running  on  a  900  MHz  Raspberry  Pi  23  can  easily  satisfy  control  loop  periods  on 
the  order  of  10_3s.  Thus,  we  assume  that  the  time  to  execute  a  piece  of  control  code,  At,  is  much 
less  that  the  control  loop  period,  T.  Figure  2  visualizes  the  parallel  execution  of  control  software 
with  system  dynamics  by  using  the  notion  of  superdense  time  [12]. 


1  https://www.raspberrypi.org/ 
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Figure  3:  An  illustration  of  how  superdense  time  models  the  execution  of  a  software  controller 
between  autopilot  sample  times.  When  a  control  loop  is  intiated  at  its  next  sample  period,  the 
code  executes  instantaneously  along  the  steps  along  the  index  set  N.  The  control  value  is  then 
written  to  the  actuator  and  held  until  the  next  sample  time. 

Superdense  time  is  the  set  S  =  R  x  N  where  R  is  a  real  valued  clock  and  N  C  N>o  is  a  finite  subset 
of  indices.  The  continuous  portion  of  system  state,  x(t),  evolves  along  the  physical  time  axis,  t. 
The  software  controller  state,  xc,  evolves  along  the  code  execution  axis,  n,  where  each  instruction 
is  a  discrete  step.  Consider  the  sample  time  kT  in  Figure  4.2,  where  a  sensor  provides  new  data  to 
the  control  software  routine.  The  computer  executes  each  line  of  the  routine  at  an  index  value  from 
the  code  execution  time  line.  When  the  control  calculation  routine  completes,  it  writes  the  actuator 
value  for  that  sample  time,  u(kT),  to  the  memory  location  used  to  set  actuator  speed  or  position. 
This  actuator  value  will  be  held  until  the  next  sample  time,  ( k  +  1  )T,  and  the  software  controller 
executes  again  along  its  discrete  axis. 

Note  that  this  model  of  software  controller  execution  is  not  limited  to  directly  interfacing 
software  and  hardware.  It  can  also  capture  higher  level  control  algorithms  that  read  (write)  from 
(to)  data  channels.  For  example,  a  multi-agent  consensus  algorithm  might  block  on  a  queue  that 
provides  new  neighbor  and  local  state  information.  The  consensus  software  would  compute  the 
resulting  control  value  as  soon  as  the  new  data  arrives  in  the  queue.  If  this  data  is  provided  at  a 
regular  interval,  the  algorithm’s  execution  would  operate  in  a  similar  way  to  the  hardware  example 
above. 


4.2  Generating  Dynamical  Models  of  Code 

The  authors  of  [22]  developed  a  graphical  model  of  software  and  created  an  algorithm  to  analyze 
it  for  runtime  behavior.  However,  there  was  no  process  to  ingest  controller  code  written  in  a  general 
purpose  language  and  automatically  create  a  model.  Figure  4.3  shows  the  process  that  translates 
source  code  written  in  C/C++  into  a  polynomial  hybrid  automata.  The  LLVM  compiler,  clang, 
compiles  the  source  le  into  LLVM’s  intermediate  representation  (IR)  bytecode.  The  intermediate 
gha  application  builds  the  polynomial  GHA  model  by  processing  a  subset  of  the  IR  instruction  set 
to  build  the  elements  in  Definition  1.  To  facilitate  analysis  of  GHA  by  future  downstream  tools, 
the  GHA  Model  output  is  generated  using  a  protocol  buffer  schema2 

2https://github.com/google/protobuf 
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IR  bytecode  les  produced  by  clang  represent  a  hierarchy  of  program  objects:  modules. 


functions,  basic  b 

Figure  4:  The  internal  process  for  TB  Generator  that  turns  controller  source  code  into  a  polynomial 
GHA. 


functions;  consequently,  each  function  has  a  list  of  instructions.  Basic  blocks  are  support 
structures  that  contain  a  sequence  of  instructions  that  have  no  branching  operations.  The  gha 
application  processes  a  single  module  by  inspecting  each  function’s  IR  instructions  in  order  of 
execution. 

We  map  sequences  of  instructions  that  do  not  branch,  or  call  other  functions,  into  a  single  mode. 
Each  program  variable  that  is  assigned  on  the  stack  becomes  a  state  variable  for  that  mode.  Using 
this  modeling  approach,  we  assume  that  there  are  no  dynamics  within  a  mode;  instead,  state 
variables  evolve  through  the  reset  functions  along  each  edge  of  the  produced  model.  The  GHA 
captures  branching  behavior  by  converting  comparison  instructions  into  guard  functions.  The 
inputs  and  outputs  of  a  software  model  are  indicated  through  function  calls  that  read  or  write 
hardware  registers.  For  example,  a  function  sync_in  ( )  might  perform  a  blocking  read  from  the 
A/D  in  Figure  4.1.  These  synchronization  functions  also  generate  the  events  contained  in  X.  All 
other  edges  in  the  model  receive  the  empty  label,  e,  indicating  that  no  synchronization  is  required 
on  that  edge. 


4.3  Code  to  Model  Example 

We  apply  the  modeling  discussion  of  the  prior  section  to  the  example  code  snippet  in  Figure  4.4, 
which  increments  an  integer  value,  a,  from  0  to  9.  A  diagram  of  the  automatically  generated  hybrid 
model  is  shown  in  Figure  4.5. 


1 


5 


int  main() { 

int  N  =  IS; 
int  a  =  0; 

for (int  i  =  9;|  i  <  N;  i++){ 
a++; 

} 

} 


Figure  4:  A  simple  for-loop  code  snippet  ingested  by  the  TB  Model  Generation  process. 


The  formal  hybrid  automata  model,  Afor,  for  this  source  code  is  generated  by  applying  Definition 
1  over  the  IR  code  that  follows  the  requirements  described  in  Section  4.2.  This  code  has  five  modes, 
Q=  {qo,---  ,<?4},  five  edges, 

E  =  {e q  =  (0,  l,  e),ei  =  (l,2,e),e2  —  (1,3,  e),  e3  =  (2,4,  e),e4  —  (4,  l,e)}  ;  and  a  state  space 
set  X  composed  of: 
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>V<?eQ 


N 


a 

i 

retval 


Note  that  the  state  variable  retval  represents  the  return  status  of  the  main( )  process.  This  model 
has  no  inputs  or  outputs,  so  U  =  Y  =  0  and  H  =  0  .  Since  there  are  no  synchronization  labels,  all 


Figure  5:  The  polynomial  hybrid  automata  that  is  automatically  generated  from  the  source  code 

in  Figure  4. 

edges  have  the  empty  label, .  Then  the  guard  function  set  is  given  by: 

d*  {0eo  <Pe1  —  N  i,  0g2  1  N,  0g3  —  1,  0g4  —  1} 

Based  on  the  assumptions  described  in  Section  4.2,  the  physical  dynamics  of  each  mode  are 
fq  =0,  e  Q.  The  software  state  variables  are  transformed  by  a  set  of  reset  functions,  R, 
associated  with  each  edge: 


Re*  — 

[10  0  0 

o]T 

Re  i  = 

[N  a  i 

retval]T 

R.^2  = 

[n  a  i 

retval]T 

R^s  = 

[N  a  +  1 

i  retval] 

R^4  = 

[n  a  i  +  1  retval] 
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Chapter  5:  Simulations  and  Abstraction 
5.1  Motivation 

The  TeamBlocks  model  generation  and  validation  process  is  illustrated  in  Figure  6.  The  TB 
Model  Generator  component  accepts  software  controller  source  code,  such  as  PID  or  consensus 
algorithms.  Using  the  LLVM4  compiler  infrastructure,  this  component  translates  C  or  C++ 
source  code  into  a  hybrid  dynamical  model,  which  is  presented  in  Section  3.1. 

The  TB  Validator  consumes  the  hybrid  automata  models  for  the  source  code  and  its  ideal 
description.  It  uses  sum-of-squares  (SOS)  optimization  [16,  20]  algorithms  to  1)  generate  an 
approximate  bisimulation  function  as  a  behavioral  certificate,  or  2)  report  that  the  code  is  invalid 
based  on  the  requirements.  Although  this  paper  is  focused  on  the  validation  of  implemented 
source  code  against  its  ideal  model,  the  TB  Validator  may  operate  over  any  pair  of  models  that 
conform  to  the  definitions  presented  in  Section  3.1. 


Figure  6:  This  figure  illustrates  the  TeamBlocks  model  generation  and  validation  process. 

The  more  precise  sense  in  which  one  system  can  be  made  behave  to  like  another  one  is  called 
approximate  (bi)simulation;  we  describe  this  next,  beginning  with  the  older  concepts  of  exact 
simulation  and  bisimulation. 

5.2  Definitions 

5.2.1  Simulation  and  Bisimulation 

Informally,  a  transition  system  73  can  simulate  a  system  T\,  if  there  exists  a  memoryless 
controller  for  73  that  makes  it  input-output  identical  to  7). In  symbols,  we  have: 

Definition  3  (Simulation).  Let  Q\,Qi  be  the  state  spaces  of  73,72,  respectively.  A  relation  S  C 
Qi  x  Qi  is  called  a  simulation  relation  of  T\  by  73  is  for  all  (q\,qi)  eS: 

•  hi(qi)  =  hiiqi) 

•  for  all  «/i  -*<r  <l[  there  exists  ®  -+<r  <?2  such  that  WiiQz)  € 


4  http://llvm.org/ 
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If  a  simulation  relation  of  72  by  T\  exists,  then  T\  is  said  to  simulate  To. 

Translating  back  to  words,  such  a  relation  S  has  the  property  that  states  can  be  associated  between 
Qi,Q2,  and  inputs  to  T2  chosen,  such  that  corresponding  outputs  of  To  are  identical  to  outputs  of  7). 
When  we  read  there  exists  in  the  above  definition,  we  can  mentally  substitute  there  exists  a 
function  to  choose,  or  there  exists  a  controller. 

Likewise,  a  bisimulation  relation  exists  between  T\  and  T2  if  T \  can  simulate  T2  and  T2  can 
simulate  Ti,  i.e.: 

Definition  4  (Bisimulation).  Let  Qi,Q2be  the  state  spaces  of  T\,T2,  respectively.  A  relation  S  C 
Q\  x  Qi  is  called  a  bisimulation  relation  between  T\  and  T2  if  for  all  (q\,qi)  e  S: 

•  hi(qi)  =  hiiqi) 

•  for  all  <7i  ~*<r  (l[  there  exists  92  ~~*o-  q'2  such  that  e 

•  for  all  <72  ^<7  <72  there  exists  <7i  ~><t  Qi  such  that  Wi-.Qz)  e  S. 

If  a  bisimulation  relation  between  T2  and  7)  exists,  then  T\  and  T2  are  said  to  be  bisimilar. 

In  other  words,  controllers  exist  to  make  either  system  input-output  identical  to  the  other. 

5.2.2  Approximate  Simulation  and  Bisimulation 

In  [8],  closely-related  definitions  were  introduced  that  allow  for  a  greater  degree  of 
approximation.  These  modified  definitions,  appropriately-enough  named  approximate 
(bi)simulation,  simply  introduce  a  distance  metric  d  and  replace  the  exact  equality  h{q\)  =  hiqo) 
by  the  approximate  equality  d(h{qi), d{q>)  )  <  e  The  definitions  follow: 

Definition  5  (Approximate  Simulation).  Let  Q\,Q2 be  the  state  spaces  of  7), 72,  respectively.  A 
relation  S  C  Q\  x  Q2  is  called  an  -simulation  relation  of  T \  by  To  is  for  all  {q\,qi)  e  S: 

•  d(h1(q1),h2(q2))  <  e 

•  for  all  q1  ->a  qx'  there  exists  q2  ->a  q2  such  that  (q^,  q2 ')  G  5. 

If  a  e-simulation  relation  of  T2  by  7\  exists,  then  is  said  to  e -simulate  T2. 


Definition  6  (Approximate  Bisimulation).  Let  Qi,Q2be  the  state  spaces  of  7), 72,  respectively. 
A  relation  S  ^  Q\  x  Q2  is  called  a  -bisimulation  relation  between  T\  and  To  if  for  all  {q\,qi)  e  S: 

•  d(h1(<q1'),h2(<q2)')  <  e 

•  for  all  q1  q x'  there  exists  q2  ~^a  q2  such  that  (qq',  q2)  G  S. 

•  for  all  q2  q2  there  exists  q1  ->a  qt'  such  that  (q^,  q2)  G  S. 

If  a  e-bisimulation  relation  between  T2  and  7\  exists,  then  Tt  and  T2  are  said  to  be  e-bisimilar. 
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As  [8]  argues,  the  importance  of  the  newer  approximate  definitions  is  that  they  enable 
additional  abstraction.  In  short,  the  relaxed  definition  may  allow  a  system  to  get  away  with  using 
fewer  states  to  approximate  a  more  complicated  system,  than  it  would  have  needed  to  exactly 
simulate  it.  For  our  purposes  in  TeamBlocks,  it  is  primarily  simulation  that  will  be  of  importance. 
We  will  look  for  controllers  that  cause  a  complicated  system  to  approximately  simulate  a  simpler 
one  and  we  will  look  for  certificates  that  the  approximate  simulation  relation  actually  exists. 

5.2.3  Simulation  Functions 

Our  certificate  that  one  system  really  does  approximately  simulate  another  one  will  take  the 
form  of  a  simulation  function.  These  are  scalar-valued  functions  of  the  states  of  the  two  systems, 
which  upper-bound  the  error  that  may  possibly  be  seen.  More  precisely: 

Definition  7  (Simulation  Function).  A  function  V :  Q\  x  Q2  is  a  simulation  function  if. 


V{qi,q2)  >  max  <*(91,92),  sup  inf  V{q[,q'2) 

V  /.  (5.1) 

If  we  have  a  function  V :  Q\  xQ2,  we  can  easily  verify  that  it  is  indeed  a  simulation  function, 
by  performing  a  local  test  at  each  joint  state  (qi,q2).  In  this  sense  a  simulation  function  is  much 
like  a  Lyapunov  function  it’s  difficult  to  find  but  easy  to  test,  and  it  implicitly  carries  global 
information  about  the  dynamical  system. 

If  (5.1)  holds  with  equality,  then  we  arrive  at  a  system  of  Bellman-like  inequalities,  and  V is 
known  as  the  directed  branching  distance.  Like  the  optimal  Value  Function,  it  can  be  computed 
by  Value  Iteration,  but  this  can  only  be  done  in  finite-state  systems,  and  the  Curse  of 
Dimensionality  makes  this  impractical  in  all  but  the  smallest  of  these.  Thus,  we  will  need  to  find 
a  better  way  to  search  for  simulation  functions.  This  will  be  the  topic  of  the  next  chapter  but 
first,  we  make  one  additional  simplification. 

5.2.4  The  Interface  Controller  and  the  Error  System 

If  we  fix  a  choice  of  controller  5  K :  Q\  x  E  Q2,  then  (5.1)  becomes, 


^(91, 92)  >  max  <*(91,92),  sup  V(q[,K(qi,a))  \ 

V  «-Wi  /.  (5.2) 

Equivalently,  we  form  a  new  system  T  by  first  (a)  forming  the  product  system  T\xT2,  then  (b) 
applying  the  feedback  interconnection  operator  with  feedback  map  K,  and  finally  (c)  defining 
the  output  function  h(q)  =  d(qi,q2),  at  which  point,  letting  Q  =  Q\  xQ2  be  the  state  space  of  T, 
(5.2)  can  be  expressed, 


5  K  is  called  the  interface  controller  because  it  causes  T2  to  satisfy  the  interface  defined  by  7i. 
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(5.3) 


V(q)  >  max  I  h(q),  sup  V{q')) 

\  q-*  <7<?' 

The  transition  system  T  is  illustrated  by  Figure  7.  In  short,  for  any  choice  of  interface  controller 
K  we  will  have  an  upper  bound  on  the  directed  branching  distance  (and  thus  a  simulation 
function),  and,  fixing  that  choice  of  controller,  the  problem  reduces  to  bounding  the  output  of 
the  error  system  T.  We  will  tackle  this  problem  with  Sum  of  Squares  optimization  in  the  next 
chapter. 


Figure  7:  The  error  system  T  is  formed  through  a  product  of  the  complicated  system  T2  with 
the  simple  (or  abstract)  system  7),  and  its  feedback  interconnection  with  an  interface 

controller  K. 
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Chapter  6:  Sum  of  Squares 

In  an  effort  to  tame  the  Curse  of  Dimensionality,  in  this  chapter  we  will  attempt  to  search  for 
polynomial  simulation  functions,  given  a  fixed  choice  of  interface  controller.  To  do  this,  we  will 
form  a  convex  relaxation  (really  a  stricter,  sufficient  problem),  that  in  principle  will  allow  us  to 
efficiently  search  for  these  polynomials  by  solving  a  convex  optimization  problem.  If  we 
succeed  in  finding  such  a  polynomial,  then  we  will  have  an  easily-checked  certificate  that  one 
system  -simulates  another.  (If  on  the  other  hand  the  optimization  is  infeasible,  then  we  will  have 
proven  nothing.)  The  method  by  which  we  search  for  these  polynomials  is  known  as  Sum  of 
Squares  programming  an  efficient  method  for  finding  positive  polynomials. 


6.1  Introduction 

A  polynomial  p  in  one  or  more  variables  is  a  sum  of  squares  (SoS)  if  it  is  of  the  form 


N 


=  Y,p? 


i=i 


(6.1) 


where  N  >  1  and  pi,..., /Mr  are  arbitrary  polynomials.  Clearly,  such  a  polynomial  is  everywhere 
positive.  What  is  more  interesting,  however,  is  that  whereas  it  is  NP-hard  to  test  whether  an 
arbitrary  multivariate  polynomial  is  pointwise  positive,  one  can  test  whether  it  is  a  sum-of- 
squares  in  polynomial  time  by  solving  a  semidefinite  program  (SDP).  Thus  SoS-ness  is,  in  the 
general  (i.e.  multivariate)  case,  a  sufficient  but  not  necessary  condition  for  positivity.  Our 
approach,  then,  will  be  to  form  systems  of  polynomial  inequalities,  to  turn  these  into  positivity 
constraints  on  polynomials,  and  to  then  to  replace  those  positivity  constraints  by  sufficient  (but 
not  necessary)  SoS  constraints.  In  fact,  the  procedure  for  converting  SoS  constraints  to  SDP 
constraints  is  relatively  simple: 


Given  an  appropriately-large  6  vector  b  =  (b\,...,bn)  of  basis  polynomials  (e.g.,  monomials  to 
sufficient  degree),  p  is  a  sum-of-squares  if  and  only  if 

p  =  bTSb  (6.2) 

for  positive-semidefinite  (PSD)  matrix  S  =  ST  >  0.  Thus,  constraints  that  polynomials  be  SoS, 
can  quickly  be  reduced  to  constraints  that  real  matrices  be  PSD. 

6.2  Simulation  Function  Search 

Our  algorithm  will  take  the  following  form: 

1.  Form  product  automaton  of  model  (in  feedback  with  interface  controller)  with 
abstraction. 


6  There  are  a  number  of  techniques  for  reducing  the  number  of  basic  elements  needed,  e.g.  by  exploiting  properties  of  the 
Newton  Polytope;  for  more  information,  see  [20]. 
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2.  Define  error  output  polynomial  for  the  product  automaton  e.g.  h2(x2)  ~  /?i(xi). 


3.  Choose  a  polynomial  basis  (e.g.,  monomials)  in  which  to  search  for  a  candidate 
simulation  function. 

4.  Generate  constraints  to  link  (the  coefficients  of)  various  polynomial-valued  decision 
variables. 

5.  Solve  a  convex  program  for  the  polynomial  coefficients. 


The  main  question  arises  at  step  4:  What  constraints  should  be  formed  to  define  a  simulation 
function  for  a  polynomial  hybrid  automaton?  And  what  decision  variables  will  we  need  to 
define?  We  answer  these  questions  next. 


6.2.1  The  Optimization  Problem 


Decision  variables:  For  each  mode  q  G  Q,  we  define  a  polynomial  decision  variable  Vq; 
together  these  describe  the  hybrid  simulation  function  (x,  q)  ■->  Vq(x).  (The  degree  of  the  Vq  is 
an  algorithm  parameter.)  Additionally,  following  the  s-procedure,  we  define  for  each  edge  e  G 
E  a  “Lagrange  multiplier”  polynomial  Se  (again  with  some  chosen  degree).  Finally,  we  also 
introduce  a  scalar  decision  variable  y  >  0  to  upper-bound  the  simulation  function. 

Mode  constraints:  For  each  mode  q  G  Q,  we  define  constraints  that  the  following  polynomials 
beSoS: 


Vq(x)  -  hq(x)Thq(pc) 


dVn 


1(x)fn(x) 


dx  v"'7<?v 
Y  ~  Vq(x) 


The  simulation  function  upper  —  bounds  the  instantaneous  error 

The  simulation  function  is  non  —  decreasing  within  a  mode 
y  upper  —  bounds  the  simulation  function 


Edge  constraints:  For  each  edge  e  =  (s,  d)  G  E,  we  define  constraints  that  the  following 
polynomials  be  SoS: 


Vs(x)-Vd(Re(x))  +  Se(x)0e(x) 


Meaning:  <£(x)  >  0  =>  V^(x)  >  Vd(Re(x)) 


Se(x) 

Optimization  Problem:  Solve, 


min  y 

{Vq}qeQi{Se}eeE'Y 

subject  to  the  preceding  constraints. 


(6.8) 


6.3  Example 

Simulation-function  certificate  generation  was  demonstrated  by  comparing  two  relatively 
simple  two-well  polynomial  systems  of  the  form, 
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X  =  V 


V  =  -—Ui(x)  -  C V  +  U 

for  polynomials  U\,Lh  (potential  functions),  scalar  C>  0  (damping  factor),  and  control  input  u 

(force).  The  system  is  illustrated  in  Figure  6.1;  in  this  case,  the  polynomials  were  given  by, 

x2(x2  -  2)  +  1 
U\{x)  = - - - 

U2(x)  =  5(x  +  l)2  -  f  (x  +  l)3  +  J^(x  +  l)4 

After  the  parallel  composition  of  Systems  1  and  2,  the  certificate  polynomial  V  was 
computed  by  sum-of-squares  optimization.  For  a  polynomial  basis,  all  monomials  of  order  less 
than  or  equal  to  8  were  used.  The  YALMIP  frontend  was  used  to  generate  SoS  constraints,  and 
the  SDPT3  was  used  to  solve  the  resulting  convex  optimization  problem. 

The  resulting  certificate  polynomial  V  of  the  joint  state  ( xi,vi,x2,v2 )  has  495  terms,  a  few  of 
which  are  shown  below,  together  with  global  upper  bound  y  on  the  error,  obtained  via  SoS 
constraints: 

F (x)  =  20.00081495  -  0.0027xi  +  0.0033vi  +  ••  •  -  \.1699e  -  05x2v27 -  4.5964e  -  06v28 
y  =  20.2882 


Figure  8:  Two  similar  polynomial  potential  wells 


A  few  things  are  worth  noting.  First,  we  are  able  to  compute  a  certificate  polynomial;  and, 
moreover,  we  obtain  a  global  upper  bound  for  the  error  between  the  two  systems  (20.2882). 
However,  note  also  how  ill-conditioned  the  problem  is  when  using  this  monomial 
representation:  Products  of 
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Figure  9:  The  simulation  function  polynomial  V  for  the  two- well  system 

g 

very  small  numbers  (e.g.  4.5964e-06)  with  large  powers  (in  this  case  vi)  appear.  This  quickly 
starts  to  present  difficulties  for  the  interior  point  solvers  (e.g.  SDPT3)  that  we  use,  and  as  a 
result  the  technique  will  require  modifications  for  use  on  larger  problems.  In  particular,  future 
research  may  investigate  the  use  of  other  polynomial  and  trigonometric  bases. 
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Chapter  7:  Modeling  Vehicles  and  Teams 

7.1  Models  Overview 

In  this  section,  we  describe  a  collection  of  continuous-time  state-space  models  that  are 
provided  as  reference  models  for  use  with  the  TeamBlocks  Framework.  The  models  provided  in 
the  current  version  are  tailored  for  TeamBlocks  Demonstration  1  (scheduled  for  December  2015) 
and  are  intended  to  form  the  basis  of  a  hierarchy  of  abstraction  in  multiagent  systems  that  spans 
low-level  vehicle  dynamics  to  higher-level  team  behaviors. 

Each  of  the  examples  that  follows  has  a  single  mode  of  operation,  which  maps  to  a  single  Mode 
in  the  TeamBlocks  schema,  without  any  edges.  The  dynamics  are  represented  by  the  Mode’s 
vector_field  field,  which  describes  a  polynomial  that  is  the  right  hand  side  of  an  ordinary 
differential  equation  (ODE)  in  the  state  and  input  variables. 

7.2  Vehicle  Models 

Here,  we  describe  vehicle  models  that  appear  in  TeamBlocks  Demonstration  1  abstraction 
hierarchy. 

We  begin  with  the  simplest,  highest-level  kinematic  descriptions  of  vehicles  as  integrators,  and 
finish  with  more  complicated  models  of  ight  dynamics,  before  proceeding  to  team  abstractions  that 
aggregate  these  models  into  larger  formations. 

7.2.1  Cartesian  Single-Integrator 

Summary  The  single-integrator  is  the  simplest  kinematic  vehicle  abstraction.  It  is  useful  in 
situations  when  dynamics  can,  relative  to  the  required  tolerances,  be  entirely  abstracted  by 
lowerlevel  controllers.  This  model  allows  the  application  of  the  extremely  well-developed  tools  of 
linear  control  theory. 

Parameters  Models  in  this  class  are  parameterized  by  the  spatial  dimension  n  e  {  23 }. 


Modes  The  system  consists  of  a  single  mode. 

Inputs  The  input  to  the  system  is  an  instantaneous  velocity  command  u  eR n=  U. 


State  Space  The  state  x  e  x  =  R"  of  the  vehicle  is  its  n-dimensional  position  in  Euclidean  space. 
Dynamics  The  position  x  e  X  evolves  according  to  the  rst-order  linear  ODE  x'  =  u. 


7.2.2  Cartesian  Double-Integrator 

Summary  The  double-integrator  is  used  to  abstract  vehicles  in  cases  when  their  acceleration 
is  bounded  but  their  orientation  is  unimportant.  As  in  the  single-integrator  case,  systems  that  can 
be  abstracted  to  this  model  can  be  treated  with  the  many  tools  of  linear  control  theory. 

Parameters  Models  in  this  class  are  parameterized  by  the  spatial  dimension  n  e  {  23 } . 
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Modes  The  system  consists  of  a  single  mode. 

Inputs  The  input  to  the  system  is  an  instantaneous  acceleration  command  u  £R n=U. 


State  Space  The  state  space  X  =  R”  x  R"  3  (p,v)  consists  of  n-dimensional  positions  and  velocities. 

Dynamics  The  state  evolves  simply  according  to  the  linear  ODE, 

p'  =  v 
v'  =  u . 


7.2.3  SE(2)  Integrator 

Summary  The  SE(2)  integrator  is  a  kinematic  model  of  a  nonholonomic  vehicle  in  two- 
dimensional  space.  It  is  used  to  abstract  ground  vehicles,  surface  vessels,  and  fixed-wing  aircraft 
constrained  to  constant-altitude  flight. 

Modes  The  system  consists  of  a  single  mode,  containing  the  ODEs  described  in  the  next 
paragraphs. 

Inputs  An  input  u  =  (v,co)  £R  2=U  to  the  system  consists  of  a  (signed)  speed  command  v,  together 
with  a  angular  velocity  command  co. 

State  Space  The  state  space  X  =  R2  xSO(2)  3  (p,R)  of  the  system  consists  of  all  2d  rigid-body 
positions  and  orientations. 

Dynamics  The  state  evolves  according  to  the  ODEs, 

p'  =  vRe  i 
R'  =  coRJ 

where 


is  the  90-degree  rotation  matrix  (named  J  for  its  analogy  to  the  imaginary  number),  and  e\  =  (1,0) 
is  the  rst  element  of  the  natural  basis  for  R2. 


7.2.4  SE(3)  Integrator 

Summary  The  SE( 3)  integrator  is  a  simple  kinematic  model  of  a  nonholonomic  vehicle  in  three- 
dimensional  space.  It  is  used  to  abstract  fixed-wing  aircraft  and  other  vehicles  (e.g.,  unmanned 
underwater  vehicles  (UUVs))  that  perform  three-dimensional  maneuvers  by  steering. 
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Modes  The  system  consists  of  a  single  mode,  containing  the  ODEs  described  in  the  next 
paragraphs. 

Inputs  An  input  u  =  (v,co)  e  RxR3  =  U  to  the  system  consists  of  a  (signed)  speed  command  v, 
together  with  a  angular  velocity  command  co. 

State  Space  The  state  space  X  =  R3  x  50(3)  ~  SE( 3)  3  (p,R)  of  the  system  consists  of  all  3d 
rigid-body  positions  and  orientations. 

Dynamics  The  state  evolves  according  to  the  ODEs, 

p'  =  vRe  i 

R'  =RQ. 

where 


o 


n  = 


wa 


— ua 

0  — Qj’I 

0 


and  e\  =  (1,0,0)  is  the  fi  rst  element  of  the  natural  basis  for  R3. 


7.2.5  Flight  Model 

Summary  This  is  a  quasi-static  model  of  6-degree-of-freedom  fixed-wing  flight  with  12  state 
dimensions.  Quasi-static  in  this  case  means  that  the  aerodynamic  forces  are  entirely  a  function  of 
the  aircraft’s  orientation  and  (angular  and  translational)  velocity,  as  opposed  to  also  depending  on 
additional  state  possessed  by  the  surrounding  fluid. 

Modes  The  flight  model  described  here  consists  of  a  single  mode,  7  containing  the  equations  of 
motion  that  are  described  next. 

Inputs  The  input  space  for  the  aircraft  is  U  =  R>o  x  R3;  an  input  u  =  ( uv,um )  consists  of, 

•  mv  e  R>o,  the  engine  thrust,  in  units  of  force,  and 

•  Uco  =  (i Ur,ue,ua )  e  R3,  the  rudder,  elevator,  and  aileron  deflections,  in  radians. 


7  We  anticipate  that  subsequent  models  will  capture  piecewise  features  of  lift  and  drag  curves  through  separate  glide  and  stall 
modes,  each  subject  to  different  aerodynamic  forces,  and  with  a  passport  (or  guard)  dependent  on  the  angle  of  attack. 
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State  Space  The  state  space  of  the  aircraft  is  X  =  R3  x  SO( 3)  x  R3  x  R3;  states  consist  of 

•  pi  e  R3,  the  inertial  position  of  the  aircraft, 

•  Rib  eSO(3),  the  rotation  from  the  aircraft’s  body  frame  to  the  inertial  frame, 

•  vb  e  R3,  the  body  velocity  of  the  aircraft,  and 

•  cob  e  R3,  the  angular  velocity  of  the  aircraft  in  the  body  frame. 

The  subscripts  I  and  B  are  used  above  to  denote  inertial  and  body  frames,  respectively. 
Additionally,  we  denote  by 


- 


0 

WB,  3 


-WB,  2 


-Wfl,  3 
0 

Wfl,  i 


WB,  2 
0 


the  cross-product  matrix  corresponding  to  cob- 

Dynamics  The  state  variables  evolve  according  to  the  ODEs, 

P]  =  Ribvb 
Rib  =  Rib^b 

vB  =  M_1(— /2BMvB  +  fv(x)  +  Gv(x)uv ) 

coB  =  +  fw  00  +  Gm(x)um) 


where, 

•  fv(x)  =  fv(.RiB>  vb>  Mb)  and  (x)  =  (R/B,  vb,  o>b)  are  aerodynamic  forces  and  moments, 

to  be  specified, 

•  M  =  ml  G  R3x3  for  m  >  0  is  the  aircraft  mass  matrix, 

•  /  =  /r  >  0  G  R3x3  is  the  body-frame  inertia  tensor  about  the  aircraft  center  of  mass 

•  Gv(x)  =  Gv(RlB,  vb,  o>b)  G  M.3x1  and  =  Gw(^/bTBiwb)  £  l3x3  are  the  decoupling 
matrices  that  relate  control  surface  deflections  to  moments. 

The  functions /v,/(o,  Gv,  and  Gw  are  required  to  be  polynomial  for  use  in  the  larger  framework,  and 
in  particular  in  the  Sum  of  Squares  optimization. 

7.3  Team  Abstractions 

When  multiple  vehicles  are  combined  into  a  team  or  platoon,  they  are  then  controlled  to  behave  as 
a  single  composite  system,  called  the  Team  Abstraction.  The  Team  Abstraction  enables 
coordinated  behavior  among  the  agents,  and,  because  it  generally  has  lower  state  dimension  than 
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the  product  of  the  vehicles’  states,  it  additionally  helps  to  abstract  or  simplify  the  state  space  for 
planning  and  verification  purposes. 


7.3.1  Formation  Abstractions 

Summary  The  formation  abstraction  treats  the  formation  as  a  single  rigid  body,  together  with  a 
number  of  shape  parameters  that  describe  degrees  of  freedom  e.g.  the  spacing  of  agents  in  a  line, 
or  the  apex  angle  of  a  V  formation.  It  is  formed  as  a  Cartesian  product  of  an  integrator  in  the  shape 
parameters,  with  an  SE(3)  integrator,  as  described  next. 8 

Parameters  A  formation  abstraction  is  parameterized  on, 

•  a  number  m  GNof  degrees  of  freedom  for  the  formation, 

•  a  set  Q  C  Rm  of  possible  shape  parameters,  and 

•  a  shape  function  g  =  (gi,...,gN)  :  Q  ==>  (R3)'v,  where  N  is  the  number  of  platforms. 

The  example  of  a  “V”  formation  is  given  later  in  this  section. 

Modes  The  system  consists  of  a  single  mode. 

Inputs  An  input  u  =  (v,co,£)  eRxR3xQ  =  (/to  the  system  consists  of  a  (signed)  speed 
command  v,  an  angular  velocity  command  co,  and  a  shape  velocity  command  c. 

State  space  The  state  of  the  formation  abstraction  is  a  tuple  x  =  (p,R,q )  e  R3xS0(3)x<2  =  X. 


Dynamics  The  formation  abstraction’s  state  evolves  according  to, 

p'  =  vRe  i 

R'  =  RQ 

q'  =  £ 


where  again 


8  More  generally,  a  formation  abstraction  can  be  formed  as  a  Cartesian  product  of  (a)  a  model  with  any  dynamics  on  the 
shape  parameters  (e.g.,  double  or  single  integrator),  and  (b)  a  vehicle  model  whose  state  space  contains  an  SE[ 3]  subspace 
(e.g.,  the  aircraft  model). 


Approved  for  Public  Release;  Distribution  Unlimited. 
24 


Outputs  The  formation  abstraction  generates  an  output^  —  +  p]{Li  €  1  —  (M3}A 


Example  The  “V”  formation  on  N  platforms,  as  e.g.  in  Figure  7.1,  is  defined  by, 
•  m  =  2  degrees  of  freedom. 


shape  parameters  W,L)  e  (Q/r]xR>o  =  Q,  that  define  the  V  angle  and  integrant  spacing, 
respectively,  and 


Figure  10:  A  formation  abstraction  for  “V”  -type  formations. 

•  the  shape  function  g  =  (gi,...,g,v) :  Q  ==>  (R3)'v,  defined  by 

9i{e,L)  =  (0,0,0) 

9k(0,L)  =  |(— cos(0/2),sin(0/2),O)  jfe  >  1  is  even 

9k(0.  L)  =  *±l(-ooB(e/2),-sin{0/2),O)  >  1  is  odd  . 

Thus  the  “V”  formation  summarizes  the  state  of  N  agents  by  a  total  of  8  state  dimensions. 
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Chapter  8:  TeamBlocks  Tools 


8.1  GHA  Schema 

The  tbautomaton  command  line  tool  interacts  with  the  TeamBlocks  Framework  through  a 
serialization  schema,  which  provides  the  concrete  realization  of  the  abstract  TB  automaton.  This 
section  describes  TeamBlocks’  first  such  schema. 

One  motivation  for  providing  a  concrete  serialization  is  that  it  allows  the  TeamBlocks 
Framework  to  support  analysis  tools  implemented  in  multiple  languages.  In  particular,  while 
TeamBlocks  now  performs  LLVM-based  code  analysis  in  C++,  its  current  version  performs 
bisimulation  function  computations  separately  in  MATLAB  using  the  SOSTOOLS  toolbox  for 
Sum  of  Squares  programming. 

Typical  choices  for  a  serialization  format  include  s-expressions,  JSON,  and  Google  Protocol 
Buffers  (protobuf).  Matlab  and  C++  both  have  industry  and  community  support  for  Google 
Protocol  Buffers  and  JSON;  by  using  one  of  these  formats,  TB  automaton  models  can  be 
imported/exported  using  mature  deserialization/serialization  libraries.  The  current  version  of 
TeamBlocks ,  including  the  the  tbautomaton  program,  therefore  now  supports  a  concrete  protobuf 
serialization  of  the  abstract  TB  automaton.  We  expect  that  other  tools  may  interact  with 
tbautomaton  by  reading  TB  automaton  objects  following  the  schema  defined  in  this  section. 

Polynomials  Bisimulation  function  computation,  using  the  Sum  of  Squares  technique,  assumes 
that  systems  and  their  controllers  are  formulated  in  terms  of  polynomial  expressions.  Thus,  we 
define  the  Expression  Schema,  shown  in  Figure  12,  together  with  its  required  components,  to 
represent  these  polynomials.  Polynomial  expressions  are  naturally  expressed  in  a  tree  structure 
where  each  node  in  the  tree  may  be  a  binary  operation,  such  as  addition  (+),  a  variable  (x),  or  a 
numeric  literal  (1.0).  For  example,  consider  the  simple  two-variable  polynomial  p(x,y )  :=  2x  +  y. 
The  tree  structure  for  this  polynomial  is  shown  in  Figure  1 1 .  The  leaf  nodes  of  the  tree  are  always 
a  variable  or  literal,  i.e.  2,  x,  and  y.  The  intermediate  nodes  represent  sub-expressions  of  the  total 
polynomial. 

Using  these  facts,  we  define  the  Expression  Schema  that  is  composed  of  one  of  the  following: 
1)  a  BinaryOperator,  2)  a  double-valued  literal,  or  3)  a  string-valued  variable.  BinaryOperators  are 
themselves  composed  of  left  and  right  operand  Expressions  as  well  as  a  mathematical  Operator: 
+,-,  or  x.  These  objects  are  all  required  to  make  a  valid  BinaryOperator  object. 


Modes  and  Edges  The  TB  automaton  definition  in  Section  3  requires  these  polynomial  Expression 
objects,  as  well  as  several  other  components.  Figure  11  shows  the  schema  for  Modes  and  Edges 
that  represent  a  TB  automaton  object.  Each  Mode  object  has  an  integer  identifier,  id,  as  well  as 


Figure  11:  A  tree  structure  representing  the  expression  2x  +  y. 
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Expression 


Binary  Ope  rat  or  operation 

double  literal 

string  variable 

Binary  Operator 


'Operator  [  ADD,  3UBr  ML'L  \ 

Operator  operator 

Expression  left_  argument 

Expression  n  gh  t_ar  g  urnem 

Figure  12:  This  diagram  shows  the  structure  of  the  Polynomial  Schema  used  in  the  TeamBlocks 
Framework.  Green  colored  blocks  are  optional  fields,  one  or  more  Expression  objects  that 
represent  the  dynamics  of  the  Mode.  Edges  maintain  their  source  and  destination  Modes  using 
the  integers  srcjd  and  dest_id  as  well  as  a  label  string  that  may  be  used  in  synchronization 
operations.  They  must  also  contain  a  passport  expression  and  one  or  more  transformation 

expressions. 

Schema  Realization  in  protobuf  format  The  current  TeamBlocks  Framework  realizes  the  schema 
concretely  using  the  Google  Protocol  Buffers  serialization  library;  the  descriptor  (.proto)  is  shown 
in  Figure  12.  The  process  that  generates  a  TeamBlocks  Automaton  in  this  format  is  described  in 
the  next  section. 

Schema  Realization  in  MATLAB  Corresponding  to  their  protobuf  representation,  GHAs  also 
have  a  native  MATLAB  struct  representation  as  shown  in  Figure  8.5. 

8.2  Functions 

8.2.1  Core  Functions 

tbproduct(tbl,  tb2) 

•  Inputs:  two  GHA  Models. 


Figure  13:  This  diagram  illustrates  the  high  level  structure  of  TB  automaton  models.  This  schema 
is  ingested  by  optimization  tools  to  perform  the  bisimulation  function  search  for  code  certification. 
Red  colored  blocks  indicate  required  fields. 


Approved  for  Public  Release:  Distribution  Unlimited. 
27 


Returns:  the  composed  GHA  model. 


•  This  function  performs  a  composition  that  synchronizes  two  GHA  on  edge  labels. 


tbsim(tb,  startModeld,  xO,  time_variable,  input,  maxTime,  maxlters,  display) 

•  Inputs:  A  GHA  model,  the  start  mode  identifier,  state  variable  initial  condition(s),  the 
time  variable.  The  input  signal,  maximum  simulation  time,  maximum  iterations,  and 
display  are  all  optional. 

•  Returns:  A  struct  of  the  simulation  history. 

•  This  function  simulates  a  GHA  model  from  initial  conditions  until  a  maximum  time  or 
number  of  iterations. 

tbloConnectftb,  output_names,  input_names) 

•  Inputs:  A  GHA  model,  output  variable  names,  and  input  variable  names. 

•  Returns:  A  new  GHA  with  the  10  variables  substituted. 

•  This  function  connects  the  outputs  of  an  GHA  model  to  its  inputs  based  on  the  given 
name  vectors.  For  example,  the  tbdemo_ghaexecution  script  uses  this  operation  to  make 
feedback  system  out  of  the  product  of  the  linear  system  and  PI  controller. 


thread)  lename) 

•  Inputs:  String  of  a  filename. 


•  Returns:  A  GHA  struct  built  from  the  protocol  buffer  binary  file. 

•  This  function  loads  a  protobin  schema  instance  and  creates  a  Matlab  representation  using 
the  structs  described  above. 
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1 


Ilf 


15 


21  f 


25 


30 


id  5 


import  'google/ p  rotobuf/descriptor.  proto'  f 
package  TbHsg? 


message  TbHsg  { 

optional  google. protobuf^FileOescripterSet  file_descripter_set  *  If 
required  node  node  ■  2; 
required  Edge  edge  ■  3f 

> 


message  Hade 
req  u  i  red 
repeated 


{ 

umtG4  id  *  If 

Expression  vector. fie Id  *  2f 


> 


message  Edge 
required 
required 
required 
req  u  i  red 
repeated 


{ 

uintfi4  src.id  •  1 
uintti4  dst.id  ■  2 
string  label  ■  3 
Expression  passport  *  4 
Expression  transfer™  *  5 


> 


message  Expression  { 

optional  BmaryOpe rator  operation  ■  1 
optional  double  literal  »  2 

optional  string  variable  *  3 

> 


message  BmaryOpe  rator  { 
enu*  Operator  { 

ACC  ■  1; 

HUL  *  2f 
SUB  »  3; 

If 

required  Operator  operator  *  i 

required  Expression  lef t_a rgui*ent  *  2 

required  Expression  right.argurient  *  3 


Figure  14:  The  Protocol  Buffer  implementation  of  the  schema  in  Figures  12  &  13. 
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%%  Mode  Struct 

mode  ■  struct( )  ? 

■ode-  id 

% 

String  identifier 

■ode- vector,  field 

% 

1 factor  of  polynomials 

node.output.f  unction 

h 

Expression/ tunc tion 

nude .  state. va  nable 

% 

Vector  of  variables 

fiade.is.accept 

h 

Boolean 

node . o  u  tgo 1 ng_ edge, idx 

% 

Integer 

%%  Edge  Struct 

edge  *  struct; S? 

edge- s  re.  id! 

% 

String 

edge.dst_id 

H 

String 

edge. Label 

% 

String 

edge.passpo  rt 

% 

Vector  of  polynomials 

edge. t ran s fora 

% 

Vector  of  polynomials 

%%  The  GHA  Struct 

gha  *  struct (J; 

g ha. node 

% 

Vector  of  modes 

gha.  edge 

% 

Vector  of  edges 

g ha. io_ variables 

Polynomials 

Figure  15:  The  MATLAB  struct  representation  of  GHAs.  tbCompileToMat(tb,  name) 


tbCompileToMat(tb,  name) 

•  Inputs:  A  GHA  struct,  and  name  for  the  output  function. 

•  This  function  takes  a  GHA  model  and  compiles  it  into  a  set  of  Matlab  functions.  One 
would  use  this  function  to  make  a  GHA  model  execution  more  performant  than  running 
tbsim  on  the  same  GHA  model. 

tbCertificate(tb,  degree,  tolerance) 

•  Inputs:  A  GHA  struct,  the  maximum  degree  of  polynomial  used  in  certification,  the 
constraint  checking  tolerance. 

•  Returns:  A  struct  with  polynomial  objects  that  represents  the  squared  norm  bound  on 
theoutput  of  the  system. 

•  This  function  computes,  if  possible,  the  polynomial  error  bound  certificate  of  the  GHA. 
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8.2.2  Model  Templates 

tbclock(rate) 

•  Inputs:  The  rate  in  Hz  for  the  clock. 

•  Returns:  A  GHA-based  clock  model. 

tbLinearSystem(  A,B,C) 

•  Inputs:  The  canonical  A,  B,  and  C  linear  system  matrices. 

•  Returns:  A  GHA-based  linear  model. 

tbIdealPi(nMeasurement,  nControl,  kp,  ki) 

•  Inputs:  The  number  of  measurement  and  control  signals,  proportional  gain,  and  integral 
gain. 

•  Returns:  A  GHA-based  model  of  the  theoretical  PI  controller. 

tbSE2Integrator(speed) 

•  Inputs:  Speed  of  the  vehicle. 

•  Returns:  A  GHA-based  model  of  an  SE(2)  vehicle. 

tbVee(N, speed) 

•  Inputs:  The  number  of  agents  and  speed  of  the  formation. 

•  Returns:  A  GHA-based  model  of  a  collection  of  SE(2)  vehicles  in  a  “vee”  formation 


8.2.3  LTL  Translation  and  Integration 

ltlT  oTb(ltl  Characterstring) 

•  Inputs:  A  string  of  characters  that  represent  an  LTL  expression. 

•  Returns:  An  executable  LTL  specification  automaton  using  the  GHA  data  structure. 
Thetbdemo_ghaexecution  script  illustrates  how  this  function  is  invoked. 
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Chapter  9:  Conclusion 


TeamBlocks  introduces  theory  and  software  tools  that  address  the  modeling,  composition,  and 
abstraction  of  hybrid  systems.  It  does  this  using  a  single,  unified  representation  the  generalized 
(polynomial)  hybrid  automaton  (GHA)  that  represents  both  continuous  systems  (i.e.,  systems 
described  by  differential  equations)  and  discrete  ones,  like  control  software.  Within  this 
framework,  one  can  compose,  interconnect,  and  execute  a  variety  of  models,  whether  generated 
automatically  from  C++  code  by  our  tool,  or  compiled  from  Linear  Temporal  Logic  (LTL) 
specifications,  or  constructed  by  hand  in  MATLAB. 

The  LLYM-based  model  generator,  in  particular,  is  a  key  product  of  the  TeamBlocks  effort. 

We  have  also  shown  that,  despite  being  restricted  to  polynomials,  GHAs  are  sufficiently 
expressive  to  describe  a  variety  of  systems  of  interest,  from  basic  models  like  linear  systems,  to 
standard  models  of  nonholonomic  vehicles,  to  actual  C++  implementations  of  controllers  and  that 
these  models  behave  as  expected.  Indeed,  the  TeamBlocks  library  currently  ships  with  a  number  of 
useful  model  templates  provided.  In  future,  relatively  straightforward  extensions  to  include 
trigonometric  factors  can  only  increase  the  expressiveness  of  GHA-type  models  (and  may  improve 
numerical  stability). 

By  compiling  Linear  Temporal  Logic  (LTL)  specifications  to  monitor  automata,  we  are  also 
able  to  perform  runtime  verification  efficiently  and  within  the  same  framework  of  automaton 
execution. 

The  most  challenging  part  of  this  work  has  been  the  development  of  automatic  proof 
techniques,  based  on  Sum  of  Squares  (SoS)  optimization,  to  find  polynomial  certificates  that  verify 
that  one  system  abstracts  another.  Here  we  have  developed  a  formulation  of  an  optimization 
problem  over  polynomials  as  a  sum-of-squares  optimization  problem,  which  in  principle  has  the 
advantage  of  convexity  and  therefore  tractability;  moreover,  we  developed  prototype  software  that 
is  able  to  search  for  these  polynomials  using  standard,  state-of-the-art  solvers.  However,  the 
resulting  optimization  problems  appear  to  be  numerically  ill-conditioned  in  practice,  which  may 
be  a  symptom  of  the  use  of  monomial  bases;  as  a  result,  the  generation  of  certificates  for  more 
complicated  systems  will  likely  require  the  use  of  other  polynomial  or  trigonometric  basis 
functions,  and  must  be  left  for  future  work. 
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